#include<iostream>
using namespace std;
#include<iomanip>
#include<fstream>
#include<cmath>

#include"library.h"

/* ISW1d:	解出1维无限深方势井本征值和波函数。假定 hbar^2/2m = a = 1，ISW1d: Infinite-Square-Well-1-dimension。解出的波函数写进"ISW1d_wave.txt"中。
 * int n:	(0,1)上的节点数
 * double *e:	本征值，从小到大排列
 */ 
double ISW1d(int n, double *e){//1d Infinite Square Well

	int i,j;
	double h = 1.0/(n+1);// 将区域 [0,1] 分为 n+1 份，所以其中有 n 个节点, 步长为h
	double *a = new double [n*n];//课件中的矩阵
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
			if(i==j){
				a[i*n+j] = 2/h/h;
			}
			else if( abs(i-j)==1 ) {
				a[i*n+j] = -1/h/h;
			}
			else a[i*n+j]=0;
		}
	}

	jacobi(n, a, e);//对角化
	sort_eig(n, a, e);//按本征值从小到大，将本征值和本征矢量排序

	// 将波函数输出到ISW1d_wave.txt
	ofstream fp("ISW1d_wave.txt");
	for(i=0;i<n+1;i++)fp<<"0.0\t"; fp<<endl;
	for(i=0;i<n;i++){
		fp<<setprecision(4)<<(i+1)*h<<"\t";//输出节点坐标：有效数字为4位
		for(j=0;j<n;j++)
			fp<<setprecision(4)<<a[j*n+i]<<"\t";
		fp<<endl;
	}
	fp<<"1.0\t";
	for(i=0;i<n;i++)fp<<"0.0\t"; fp<<endl;
	fp.close();

	delete [] a;//释放动态内存

	return 0;
}

int main(){

	int n=5, i, j;
	double y, dy;

	double x_data[3], eig[10][3];//用作Nevill外推的输入量

	for(i=0;i<3;i++){
		x_data[i] = 1.0/(n+1);//步长
		double *e = new double [n];
		ISW1d(n, e);
		cout<<"e:"; for(j=0;j<10;j++)cout<<e[j]/M_PI/M_PI<<","; cout<<endl;
		for(j=0;j<10;j++) eig[j][i] = e[j];//本征值
		delete [] e;
		n*=2;
	}

	cout<<"Extrapolation:"<<endl;
	for(i=0;i<10;i++){
		Nevill(3, x_data, eig[i], 0.0, y, dy);// y为外推得到的 E_{i+1}
		cout<<"E"<<(i+1)<<"/pi/pi = "<<setiosflags(ios::fixed)<<setprecision(3)<<y/M_PI/M_PI<<endl;
	}
	return 0;
}
